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Abstract 

We  study  estimators  for  the  variance  parameter  of  a  stationary  process. 
The  estimators  are  based  on  weighted  Crauner-von  Mises  statistics,  and  certain 
weightings  yield  estimators  that  are  “first-order  unbiased”  for  o^.  We  derive  an 
expression  for  the  asymptotic  variance  of  the  new  estimators;  this  expression  is 
then  used  to  obtain  the  first-order  unbiased  estimator  having  the  smallest  variance 
among  fixed-degree  polynomial  weighting  functions.  Although  our  work  is  based 
on  asymptotic  theory,  we  present  exact  and  empirical  examples  to  demonstrate  the 
new  estimators’  small-sample  robustness. 
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1  Introduction 

Consider  a  stationary  process  Fi,  >3, . . . , with  mean  fi.  Such  processes  are  often  en¬ 
countered  in  the  context  of  steady-state  simulation.  For  instance,  the  F’s  might  represent 
successive  customer  transit  times  in  a  complicated  queueing  system  that  has  been  run  to 
steady  state.  If  one  is  interested  in  estimating  n,  the  obvious  unbiased  point  estimator 
is  the  sample  mean  A  measure  of  the  sample  mean’s  precision  is  Var(y„),  which  is 
unknown. 


1 


In  this  article,  we  investigate  new  estimators  for  Var(F„),  or  equivalently,  for 
=  n\/ar(y„).  A  related  quantity  is  also  of  interest — the  variance  parameter, 
<T^  =  limn->oo<7n.  The  literature  studies  many  variance  estimation  methods;  e.g.,  batch 
means,  independent  replications,  spectral  analysis,  overlapping  batch  means,  regenera¬ 
tion,  autoregressive  modeling,  and  standardized  time  series  (STS)  (see  [3]).  The  esti¬ 
mators  presented  herein  are  based  on  weighted  functionals  of  standardized  time  series. 
We  shall  show  that  the  new  estimators  are  asymptotically  unbiased  for  and  that  they 
have  lower  variance  than  competing  estimators. 

We  first  give  some  background.  The  standardized  time  series  is  defined  as 


Tn(t)  =  for  0  <  f  <  1, 

<Ty/n 

where  Yj  =  Y!,k=i^kljy  j  =  •••,«*  and  [J  is  the  greatest  integer  function 

(Schruben  [20]).  Under  mild  conditions  (see  Foley  and  Goldsman  [8],  Glynn  and  Iglehart 
[9],  or  Schruben  [20]),  one  can  show  that 

(v/;r(F„-/i),arn)  =>  (aW(l),(rB), 

where  W  is  a  standard  Brownian  motion  process,  B  is  a  standard  Brownian  bridge 
process  on  [0, 1],  and  =>  denotes  weak  convergence  (as  n  becomes  large)  on  D[0, 1],  the 
space  of  right-continuous  functions  on  [0,1]  having  left-hand  limits.  It  is  well-known 
that  ail  finite-dimensional  joint  distributions  of  B  are  normal  and  Cov(B(s),5(f))  = 
min(s,t)(l  -  max(s,t)),  0  <  <  1.  Further,  one  can  express  B{t)  =  tW(l)  —  W(t);  so 

it  is  easy  to  show  that  W(l)  and  B  are  independent,  and  thus  v/n(^n  ~  p)  and  aTn  are 
asymptotically  independent. 

The  remainder  of  the  paper  is  organized  as  follows.  §2  reviews  the  STS  weighted  area 
estimator  for  <7^;  the  weighted  area  estimator  will  serve  as  a  benchmark  for  comparison 
in  the  subsequent  sections.  §3  presents  new  estimators  similar  to  weighted  Cramer-von 
Mises  (CvM)  statistics,  and  establishes  some  of  their  properties.  In  particular,  we  find  a 
class  of  CvM  estimators  that  is  ‘‘first-order  unbiased”  for  these  estimators  also  have 
lower  variance  than  that  of  the  weighted  area  estimator.  Perform2mce  of  the  estimators 
is  studied  in  §§4  and  5,  where  we  present  exact  and  empirical  results,  respectively.  §6 
summarizes  and  discusses  the  results  of  the  previous  sections.  §7  proposes  a  number  of 
augmentations  to  the  basic  CvM  estimator  and  provides  conclusions. 


2  The  Weighted  Area  Estimator 


We  start  with  a  discussion  of  the  so-called  weighted  area  estimator  for  <7^,  first  popularized 
by  Schruben  [20].  Suppose  we  define 

AU:n)  . 

n 

and 

Aif)  =  [' m<rB{t)dU 

Jo 

where  /(<)  is  continuous  on  [0, 1],  not  dependent  on  n,  not  identically  zero,  and  normalized 
so  that  Var(y4(/))  =  <r^,  i.e., 

VarQ‘/(()e(l)<ii)  =  2^  £f{‘)ms{l-t)didt  =  1.  (I) 

(The  above  expression  can  be  simplified  a  great  deal;  see  [12].)  Thus,  A{f)  ~  <TNor(0, 1). 
Goldsman  and  Schruben  [13]  (also  see  [7],  [9],  and  [20])  show  that  A{f\n)  ^  A(f), 
where  denotes  convergence  in  distribution  as  n  becomes  large.  So  by  the  continuous 
mapping  theorem  (see  Theorem  5.1  of  Billingsley  [2]), 

A»(/;n)  i  A\f)  ~  c\l 


Since  the  limiting  random  variable  A{f)  is  the  weighted  area  under  a  Brownian  bridge 
process,  we  refer  to  i4’(/;  n)  as  the  weighted  area  estimator  for  <7*. 

Before  stating  the  main  result  on  the  weighted  area  estimator,  we  introduce  notation 
that  will  be  useful  in  the  sequel.  Define  the  covariance  function  Rk  =  Cov(yi,  n+k), 
k  =  0, ±1,±2, ...,  and  the  quantities  7  =  ^  ^  = 

Jo  /o/(5)dadt.  Further,  the  notation  p(n)  =  0(9(n))  means  that  |p(n)/q(n)|  <  K  as 
n  — *  00  for  some  constant  K,  and  p(n)  =  o{q{n))  means  that  p{n)/q{n)  — f  0  as  n  — »  00. 


The  main  theorem  for  the  weighted  area  estimator  is  as  follows. 

Theorem  1  (see  [8]  and  [12])  Under  mild  conditions  on  the  covariance  function,  the 
expected  value  of  the  weighted  area  estimator  is 


E[A’(/:n)|  =  <r’  +  +  o(l/n)- 


Under  an  additional  uniform  integrability  assumption  (see  Billingsley  [2]’s  Theorem  5.4 
and  its  preceding  comments),  the  asymptotic  variance  of  the  weighted  area  estimator  is 
Var(>42(/))  =  Var((7>x?)  =  2*7^ 


A .  ?.  .1  Ijty 


£v-»ii  anO/or 

3l8t 

ft'* 

IP _ 

□  □ 


Example  1  The  expected  value  of  the  weighted  area  estimator  with  constant  weighting 
function  fo{i)  =  >/l2,  for  all  t  €  (0, 1],  is  E[y4^(/o;  n))  =  <7*  +  37/n  +  o(l/n).  □ 

Henceforth,  if  the  bias  of  an  estimator  for  some  parameter  is  o(l/n),  we  shall  say  that 
the  estimator  is  first-order  unbiased  for  that  parameter. 

Example  2  If  the  weighting  function  f{t)  satisfies  F  =  7*  =  0  (in  addition  to  the 
normalizing  condition  (1)),  then  Theorem  1  says  that  E[y4*(/;n)]  =  +  o(l/n),  i.e., 

A\f;n)  is  first-order  unbiased  for  Examples  of  weighting  functions  yielding  first-order 
unbiased  estimators  for  are  /2(t)  =  \/§40(3f^  —  3t  +  1/2)  and  f{t)  =  v^)ricos(27ri<), 
*  =  1,2,...  (see  [8]).  □ 

3  The  Weighted  Cramer- von  Mises  Estimator 

In  this  section,  we  propose  several  estimators  for  based  on  different  Brownian  bridge 
functionals.  To  parallel  the  discussion  of  §2,  we  define 


and  ^ 

C(j)  s  /  gimaB(t))Ut. 

JO 

where  g{t)  is  a  weighting  function  normalized  so  that  E[C(^)]  *  <r*. 
In  the  sequel,  we  will  require  a  number  of  assumptions  to  bold. 


Assumptions 

1.  The  process  Y\,Y2,...  is  stationary. 

2.  The  constants  /x  and  <7*  satisfy  =>  oW,  where 


^n(0 


lntj(TLn«j  -m) 


3.  E£-oo  ii*  =  >  0. 

4.  ESliik’|/?*|<oo. 

5.  g"{t)  is  continuous  and  bounded  on  [0,1]. 

6.  /o  ^(0^(1  -  t)di  =  1  (normalizing  assumption). 
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Remark  1  Assumptions  1-4  are  conditions  on  the  underlying  stochastic  processes. 
Glynn  and  Iglehart  [9]  list  various  sets  of  sufficient  conditions  for  Assumption  2  to  hold; 
these  conditions  usually  involve  moment  and  mixing  conditions.  Assumptions  3  and  4 
hold  for  a  wide  variety  of  stochastic  processes.  Assumptions  5  and  6  are  simply  conditions 
on  the  weighting  function.  □ 

Under  the  Assumptions,  the  continuous  mapping  theorem  implies  that  C{g\n)  ^ 
C{g).  Notice  that  the  limiting  functional  C{g)  is  the  weighted  area  under  the  square  of 
a  Brownian  bridge;  by  way  of  contrast,  the  weighted  area  estimator’s  limiting  functional 
A^if)  is  the  square  of  the  weighted  area  under  a  Brownian  bridge. 

The  distribution  of  C{g)  with  constant  weighting  function  po(0  =  for  *'■^1  ^  €  [0> 
wais  given  by  Anderson  and  Darling  [1)  and  Smirnov  (22).  Over  sixty  years  ago,  Cramer 
[4]  and  von  Mises  [23]  studied  statistics  nearly  of  the  form  of  C{go;  n)  for  the  special  case 
of  independent  and  identically  distributed  (i.i.d.)  —  Anderson  and  Darling  [1] 

examined  the  distribution  of  C{g)  with  weighting  function  gAD{t)  =  [t(l  -  f)]"*  (which 
does  not  quite  meet  our  continuity  assumption).  However,  the  distribution  of  C{g)  with 
an  arbitrary  weighting  function  has  not  been  explicitly  determined;  see  Durbin  [6]  for 
additional  details.  With  this  historical  perspective  in  mind,  we  call  C{g]n)  the  weighted 
Cramer-von  Mises  (CvM)  estimator  for  <7^. 

If  we  observe  that 

C(a;n)  =  L^±a{^)k\7l-27,7,  +  yl). 

then  we  can  give  an  easy  0(n)  algorithm  to  calculate  C(p;  n): 

Z.,Si,S2  *—  0 
FOR  it  =  1  TO  n 
Z^Z  +  Yk 
Si  <- 5, +s(i)*Z 
S, -S,  +  j(i)Z2 
Z^Z/n 

C(a:  n)  ^  -  22S.  +  S, 

For  now,  we  will  be  interested  in  moments  of  the  CvM  estimator.  Our  main  theorem 
expresses  the  expected  value  of  the  CvM  estimator  C(g;n)  in  terms  of  its  weighting 
function  g(t)  and  the  covariance  function  R*.  In  what  follows,  we  define  G  =  fo  g(t)di. 

Theorem  2  Under  Assumptions  1  through  6, 

ElC(s;n)l  =  (T>  +  2(G-l)  +  <Kl/n). 

n 
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Proof  See  the  appendix.  □ 


Assumptions  3  and  4  allow  us  to  derive  a  useful  relationship  between  and  a\  (cf. 
Schmeiser  and  Song  [19]): 


A  corollary  of  the  main  theorem  that  gives  the  analogous  expression  for  the  bias  of 
C{g',n)  as  an  estimator  for  follows  immediately  from  Equation  (2). 

Corollary  1  Under  Assumptions  1  through  6, 

E[C(j;n)l  =  <t;  +  1(G  -  2)  +  o(l/n), 
n 

Some  examples  illustrate  the  consequences  of  Theorem  2  and  Corollary  1.  The  sim¬ 
plest  example  assumes  a  constant  weighting  function. 

Example  3  Theorem  2  implies  that  the  CvM  estimator  with  constant  weighting  function 
go{t)  =  6  has  expected  value  £((7(^0;  w)]  =  <r^ ■i-5n/n  +  o{l/n)  =  a* -f  47/n-i-o{l/n).  □ 

If  G  =  1  (subject  to  the  constraints  of  Assumptions  5  and  6),  Theorem  2  implies  that 
the  bias  of  C(g;n)  as  an  estimator  of  is  (^l/n).  In  this  case,  C(g\n)  is  first-order 
unbiased  for  <r^.  Indeed,  it  is  possible  to  give  such  a  weighting. 

Example  4  Consider  the  quadratic  weighting  function  ^2;c(0  =  51  -  c/2  +  c<  -  150t’, 
where  f  €  [0, 1]  and  c  is  real.  Theorem  2  implies  that  E[C(i?2;c;  n)]  =  -I-  o(l/n).  □ 

Similarly,  if  G  =  2  (subject  to  the  constrsunts  of  Assumptions  5  and  6),  Corollary  1 
implies  that  C{g\  n)  is  first-order  unbiased  for  <r^. 


n\/ar(F„) 

lt±Coy(Y.X) 

"  isl  jsl 


n-1 


Ro  +  2Z(l--)^ 

i=l  ” 

^  o  n 

tsi 

00  9  00  00  • 


ial 


<7*  +  -  +  o(l/n). 
n 


(2) 
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Example  5  Consider  the  quadratic  weighting  function  g2x{t)  =  42  —  c/2  +  ct  —  I20t’, 
where  t  €  (0,1)  and  c  is  real.  Since  G  =  Jog2-.c(^)dt  =  2,  Corollary  1  implies  that 
E(C(i2;c;n)]  =  +  o(l/n).  □ 

The  choice  of  weighting  function  g{t)  clearly  affects  the  variances  of  C{g,n)  and 
C{g).  (The  choice  of  weighting  function  f{t)  affects  the  variance  of  but  it  does 

not  affect  the  variance  of  the  area  estimator's  limiting  functional,  in  which  case 

\/ar(A^(/))  =  2a*.)  To  see  how,  we  shall  calculate  Var(C(ij)).  First,  we  have 

Lemma  1  (Patel  and  Read  [16,  p.  309})-  If  R  and  5  are  jointly  normal,  then 
Cov(R^5»)  =  2Cov’(/^,S). 

This  immediately  yields  the  following  theorem  on  the  limiting  variance  of  C{g\n). 

Theorem  3  In  addition  to  Assumptions  1  through  6,  suppose  that  the  C^(y;n)’s  (n  = 
1,2,...)  are  uniformly  integrable.  Then  Var(C(g;n))  —*  Var(C(g)),  where 

Var(C(5))  =  a*  T  T  g{s)9{t)CoM{B\s),B\t))dsdt 
Jo  Jo 

=  2a*  I  f  g{3)g(t)Coy^{B{s),B{t))dsdt 
Jo  Jo 

(by  Lemma  1) 

=  4a*  f  git){l  —  t)^  f  g{3)s^dsdi. 

Jo  Jo 

Example  6  Consider  the  constant  weighting  function  go(t)  6  from  Exzunple  3.  The¬ 
orem  3  implies  that  Var(C(5o))  =  4a* jb.  This  limiting  variance  is  significantly  smaller 
than  that  of  the  area  estimator,  for  which  Var(A(/))  =  2<7^  (Theorem  1).  □ 

Example  7  Consider  the  quadratic  weighting  function  ^3;c(0  Example  4;  this 
weighting  function  yields  a  first-order  unbiased  estimator  for  a^.  Theorem  3  gives 
Var(C(si2;c))  =  (c^  — 300c+26856)<T^/2520,  a  quantity  that  is  minimized  by  the  weighting 
function  g\{t)  =  52!i5o(0>  whence  Var(C(^J))  =:  121<t^/70.  This  limiting  variance  is  larger 
than  that  of  the  CvM  estimator  using  the  constant  weighting  function  go{t)  (Example  6); 
of  course,  the  estimator  for  a^  based  on  go{t)  is  somewhat  biased  (Example  3).  □ 

Example  8  For  completeness,  consider  the  quadratic  weighting  function  y2;c(0  from 
Example  5;  this  weighting  function  yields  a  first-order  unbiased  estimator  for  a'^.  The¬ 
orem  3  gives  Var(C(52;e))  =  (c^  —  240c  -f-  18144)<r^/2520.  This  variance  is  minimized  by 
the  weighting  function  ^2(0  =  ^j;i2o(0>  which  case  Var(C(pJ))  =  b2a*/35.  Although 
this  limiting  variance  is  lower  than  that  of  the  CvM  estimator  using  ^2(0i  minimum- 
variance  first-order  unbiased  quadratic  weighting  function  (Example  7),  it  is  still  larger 
than  that  of  the  unweighted  estimator  using  ^o(0  (Example  6).  □ 
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Ideally,  we  would  like  to  choose  a  weighting  function  that  minimizes  the  variance  of 
the  CvM  estimator  for  <t^  while  satisfying  the  first-order  unbiasedness  and  normalizing 
constraints;  i.c  Bnd  g(t)  that  minimizes  Var(C(j;))  subject  to 


G  =  1  = 

Jo 


With  this  goal  in  mind,  suppose  that  g{t)  can  be  written  as  an  m-degree  polynomial  in 
t,  i.e., 

5„(t)  =  t€  [0,1], 

isO 

for  coefficients  c^,  Ci , . . . ,  Cm  &nd  fixed  m.  After  some  algebra,  the  problem  becomes  that 
of  finding  the  coefficients  that  minimize 


subject  to 


m  ^  jn 

^  ^ —  =  1  and  ^ -  =  1. 


We  can  use  Lagrangian  multipliers  to  solve  the  above  system.  Here  the  Lagrangian 
is  given  by 

£(Co,  Cl ,  .  .  .  ,  Cmi  ,  ^j) 

where  Ai  and  A]  are  constants.  One  takes  the  m  -f  3  partial  derivatives  of  £,  sets  the 
resulting  equations  to  zero,  and  solves  the  resulting  system  of  linear  equations  for  the 
m  -)-  3  unknown  coefficients. 

Example  9  It  is  easy  to  show  via  the  Lagrangian  method  that  the  optimal-variance, 
first-order  unbiased,  quadratic  and  cubic  polynomial  weighting  function  is  y|(t)  =  -24-1- 
150t  —  150t^,  the  choice  studied  in  Example  7.  The  best  quartic  turns  out  to  be 

-1310  19270f  25230f’  16120t3  8060f^ 

»,«)  s  —+-5; - r-+-i - 3“' 

in  which  case  Var(C(^4))  =  1.042<7^.  We  can  go  further.  For  example,  the  polynomial 
weighting  function  of  degree  m  =  6  that  minimizes  Var(C(j7e))  subject  to  the  constraints 
(3)  is  given  by  gl{t)  =  where  the  Cj’s  are  as  follows. 
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This  choice  of  weights  yields  the  optimal  \/ar(C(ir|))  =  0.8093<t^,  which  is  comparable  to 
the  variance  of  the  unweighted  (albeit  biased)  estimator  C{go,n).  □ 

Remark  2  In  order  to  achieve  further  variance  savings,  we  can  continue  to  increase  the 
degree  of  the  polynomial  weighting  function.  However,  the  magnitudes  of  the  resulting 
coefficients  become  quite  large,  and  one  must  be  careful  to  avoid  round-off  error  as  well 
as  deleterious  second-order  effects  for  small  sample  sizes.  □ 

4  Some  Analytical  Examples 

This  section  presents  exact  analytical  results  involving  specific  stochastic  processes.  We 
shall  first  obtain  some  useful  expressions  for  the  expected  values  and  variances  of  the 
area  and  CvM  estimators.  We  assume  in  the  sequel  that  Assumptions  1  through  6  are 
still  in  effect. 

We  begin  with  an  intermediate  result  on  the  area  estimator. 

EM^(/;«)]  =  Var(A(/;n)) 

n  n  n  n  n 

=  ^i;i;/(^)/(^)vXoy(F,-7„F,-F,).  (4) 

Further,  if  A(/;n)  is  normal,  then  Lemma  1  implies  that 

V»r(A’(/;n))  =  2Var=(/l(/;n))  =  2(ElA>(/;n)l)'.  (5) 

The  analogous  result  on  the  expected  value  of  the  CvM  estimator  is  derived  next. 

n  n  n 

[p  ip 

=  TSs'<^v.r(r.(;)) 

"  Si  "  " 

=  4i!#(-)*’''«(F.-f.).  (6) 
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In  addition  to  the  standing  Assumptions,  suppose  that  Yi,Y2, . . .  ,Yn  are  jointly  nor¬ 
mal.  Then 

V.r(C(j;n))  = 

n*  n  n  n  n 

9—4  n— 1  n-1  •  •  •  • 

(by  Lemma  1) 

=  A££^(^)i7(^)*VCov^(^n -F„y„ -Fj).  (7) 

We  now  have  at  our  disposal  the  machinery  to  study  specific  examples  in  which  we 
calculate  the  exact  expected  values  and  variances  of  A^(f\n)  and  C(^;n)  for  various 
weighting  functions. 

In  particular,  for  the  remainder  of  this  section,  we  shall  work  with  a  first-order  moving 
average  [MA(1)]  process,  F+i  =  Bci  +  Cj+i,  t  =  1,2, . . .,  where  the  c^’s  are  i.i.d.  Nor(0, 1); 
so  /2o  =  1  +  R±i  =  6,  and  Rk  =  0,  otherwise.  One  can  derive 


and 


j2  j  y 

,2 


Cov(F„F*)  =  J  +  ^  for  fc<;, 


(8) 


(9) 


where  =:  (I  -1-^)*  and  7  =  —26.  Note  that  for  the  MA(1)  process,  Equation  (8)  implies 
that  =  <7^  4-  J,  a  result  that  agrees  with  Equation  (2). 

Example  10  We  concentrate  here  on  some  area  estimator  expectation  results  for  the 
MA(1)  process.  For  the  constant  weighting  function  /©(f)  =  \A2,  Equations  (4),  (8), 
and  (9)  show  that 

EIAU;")I  =  +  +  ^ +  ‘>(1/’'). 

n  n*  n 

as  implied  by  Example  1.  For  the  first>order  unbiased  weighting  function  /3(f)  = 
v'^840(3f’  ~  3/  +  1/2),  we  have  (also  see  |8]) 

<r^(2n*  +  7n*  4-  63n»  -  72)  4-  2l7(2n^  -  5n^  4-  lOn^  4-  5n  -  12) 

2n® 

7(<r’  +  67)  . 


”  2n2 

=  -j-oilfn). 


+  0{n-^) 


which  is  indeed  first-order  unbiased  for  and  thus  is  in  accord  with  Example  2.  □ 


10 


Example  11  We  consider  CvM  expectation  results  for  the  MA(1)  process.  For  the 
constant  weighting  function  go{t)  =  6,  Equations  (6),  (8),  and  (9)  show  that 


_  -2  .  S-y  +6^  .  7 

c[C(iyoin)l  —  + - 5 —  +  —5 


n-' 


=  <7^+ — +  o(l/n), 
n 

as  implied  by  Example  3.  For  the  quadratic  weighting  function  g2-A^)  =  51  —  c/2  +  cf  — 
150t*,  we  get 

cr^/  -  _  M  +  4n2  -  5)  .  ')r(24n=*  -  29n2  +  5) 

tl^^W2;c,  n)j  =  - - -  + 


=  ,>  +  l(£!+il)+0(n-) 


71° 


=  <r^  +  o(l/n). 


This  demonstrates  that  g^A^)  yields  a  first-order  unbiased  estimator  for  <7*  (independent 
of  the  choice  of  c),  as  implied  by  Example  4.  For  the  optimal  quartic  weighting  function 
gA^)^  we  have  (after  a  great  deal  of  algebra) 

655(<t^  -1-  67) 


E(C'(<7i;n)]  =  <7^-1- 


63n2 


-l-0(n"^)  =  <7^ -I- 0(1 /n), 


which  shows  that  this  weighting  function  yields  a  first-order  unbiased  estimator  for  <7*,  as 
anticipated  by  Example  9.  Similarly,  the  optimal  sixth-degree  weighting  function  g^O 
gives 

E(C'(96i”)l  =  + +  =  <7*  +  o(l/n), 

n* 

so  that  this  weighting  function  also  produces  a  first-order  unbiased  estimator  for  <t^. 
Notice  that  the  n~^  term  in  the  expression  for  E[C(pg;n)]  is  quite  large  compared  to 
the  terms  from  the  estimators  incorporating  quadratic  or  quartic  weights;  in  fact,  as 
alluded  to  by  Remark  2,  large  values  of  n  are  required  before  the  second-order  term  in 
E[C(y|;n)]  becomes  insignificant. 

For  completeness,  we  give  results  on  the  quadratic  weighting  function  ^2;c(0  =  42  — 
c/2  -f  ct  —  120t^.  In  this  case,  we  find  that  (after  some  algebra) 


E[C(^2,c;n)]  = 


<T*(n^  -t-  3n*  —  4)  7(n^  -h  18n®  —  23n’  -I-  4) 


n’ 


=  <,.+2  +  2(fl+M  +  o(„-3) 

n  n* 

=  <^n  +  o(l/n). 


which  demonstrates  that  this  weighting  function  produces  a  first-order  unbiased  estimator 
for  (independent  of  the  choice  of  c),  as  implied  by  Example  5.  □ 
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Example  12  Since  the  MA(1)  is  a  jointly  normal  process,  we  can  easily  derive  area 
estimator  variance  results  for  it.  Equation  (5)  and  Example  10  imply  that  for  weighting 
functions  fo{t)  =  VT2  and  /2(f)  =  \/840(3t^  —  3t  +  1/2),  we  have  Var(/4*(/;n))  = 
'2a*  +  0(1).  These  results  make  sense  in  light  of  Theorem  1.  □ 

Example  13  We  examine  the  variance  of  the  CvM  estimator  for  various  weighting  func¬ 
tions.  For  the  constant  weighting  function  ^o(0  =  Equation  (7)  gives  us  (after  tedious 
but  straightforward  algebra) 

Var(C(s„in))  =  X  + 

as  implied  by  Example  6.  For  the  (variance-optimal  and  hrst-order  unbiased  for  <t^) 
quadratic  weighting  function  g\{t)  =  —24  -h  1501  —  1501*,  some  algebra  yields 

Var(C(5:;n))  =  1.7286<t" -I- +  0(n-*)  =  1.7286<7" -I- o(l), 

n 

as  implied  by  Examples  7  and  9.  For  the  (variance-optimal  and  Brst-order  unbiased  for 
<7*)  quartic  weighting  function  $4(1),  we  can  obtain 

Var(C(i?*;n))  =  i.0418<tV  +  0(n-’)  =  1.0418(7^  +  o(l), 

n 

as  implied  by  Example  9. 

Finally,  for  the  (variance-optimal  and  first-order  unbiased  for  a*)  quadratic  weighting 
function  5j(l)  =  -18  -h  1201  -  1201*,  some  algebra  yields 

Var(C(5;;n))  =  1.4857(7''  +  +  0(n-*)  =  1.4857(7*  +  o(l), 

n 

as  implied  by  Example  8.  □ 

We  see  from  the  above  examples  that  the  area  and  CvM  estimators  behave  as  adver¬ 
tised  on  the  simple  analytical  MA(1)  example.  The  CvM  estimator  using  the  Anderson- 
Darling  weighting  function  (which  fails  to  satisfy  some  of  the  Assumptions)  does  not 
behave  so  nicely. 

Example  14  To  complete  our  series  of  examples  with  the  MA(1)  process,  we  consider 
the  expected  value  of  the  Anderson-Darling  estimator,  i.e.,  the  CvM  estimator  with 
weighting  function  gAoit)  =  (<(1  —  0]”*-  Then  it  can  be  shown  that 

E(C(,^in)|  =  >) 

«  a*(l  -  i)  _  2(1  -  i  -  2(/n(n  -  1)  +  c,)) 
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where  c*  fa  0.577216  is  Euler’s  constant.  Although  this  estimator  is  asymptotically 
unbiased  for  o*,  the  convergence  rate  of  the  expectation  to  is  much  slower  than  those 
of  the  previous  examples.  □ 

We  resort  to  Monte  Carlo  simulation  in  the  next  section  to  empirically  evaluate  the 
performance  characteristics  of  the  various  estimators  on  more  complicated  stochastic 
processes. 

5  Empirical  Examples 

In  this  section,  we  present  empirical  examples  illustrating  the  performance  characteristics 
of  the  following  variance  estimators; 

•  A*(/o;  n)  —  unweighted  area  estimator. 

;  n)  —  first-order  unbiased  quadratic  area  estimator  for 

•  C{go;  n)  —  unweighted  CvM  estimator. 

•  (^{9AD\n)  —  Anderson-Darling  estimator. 

•  ”)  —  minimum- variance  first-order  unbiased  quadratic  CvM  estimator  for  <r*. 

•  C{gl\n)  —  minimum- variance  first-order  unbiased  quartic  CvM  estimator  for  a*. 

•  C'(56;”)  —  minimum- variance  first-order  unbiased  sixth-degree  CvM  estimator  for 

(7^ 

•  C'(55;n)  —  minimum- variance  first-order  unbiased  quadratic  CvM  estimator  for 

These  examples  involve  the  Monte  Carlo  simulation  of  a  number  of  stationary  stochas¬ 
tic  processes: 

•  The  first-order  autoregressive  process  [AR(1)],  >^>1  =  t  =  1,2, ...,  where 

the  Ci’s  are  i.i.d.  Nor(0, 1  —  with  — 1  <  ^  <  1. 

•  The  first-order  exponential  autoregressive  process  [EAR(l)],  V'i.fi  =  <f>Yi  -I-  e,>i, 
t  =  1,2,...,  where  the  c,’s  are  i.i.d.  exponential(l)  with  probability  1  —  ^  and  0 
otherwise,  and  where  0  <  ^  <  1.  (See  Lewis  [14]  for  more  details.) 

•  The  M/M/1  queueing  system’s  waiting- time  process. 
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Tab]e  1:  Estimated  Expected  Values  of  Various  Variance  Estimators  —  AR(1),  ^  =  0.9. 
(Note  that  <t^  =  19.0  for  this  process.) 


n 

A^(,/o;n) 

AHh,n) 

Cigo;  n) 

CigAD.f*) 

C(gi,n) 

C(9i,n) 

C(p$;n) 

- -5 - 

il 

IBI 

IMl 

|K«Fn 

HueSOI 

0.316 

(0.001) 

0.294 

(0.001) 

0.222 

(0.001) 

0.298 

(0.001) 

3.5245 

8 

im 

IHI 

Bin 

IHI 

■cBI 

Bin 

0.892 

!<'■  "'04) 

6.1855 

16 

IHI 

■(iloul 

Bitl 

9.8346 

32 

1 

4.934 

(0.019) 

Kwnn 

Bin 

bbh 

mm 

64 

nil 

7.812 

(0.026) 

11.925 

(0.050) 

9.903 

(0.034) 

8.645 

(0.028) 

11.377 

(0.046) 

16.1908 

128 

14.914 

(0.067) 

13.148 

(0.044) 

11.476 

(0.035) 

16.119 

(0.066) 

14.395 

(0.047) 

12.811 

(0.039) 

15.525 

(0.062) 

17.5938 

256 

16.836 

(0.076) 

18.017 

(0.081) 

15.752 

(0.049) 

14.236 

(0.040) 

17.968 

(0.074) 

17.207 

(0.055) 

16.073 

(0.046) 

17.525 

(0.068) 

18.2969 

512 

17.989 

(0.081) 

18.824 

(0.085) 

17.349 

(0.052) 

16.164 

(0.043) 

18.783 

(0.078) 

18.507 

(0.059) 

17.924 

(0.051) 

18.496 

(0.072) 

18.6484 

1024 

18.084 

(0.053) 

17.285 

(0.044) 

18.838 

(0.078) 

18.867 

(0.060) 

18.617 

(0.052) 

18.687 

(0.072) 

18.8242 

2048 

18.799 

(0.084) 

18.599 

(0.053) 

18.073 

(0.045) 

19.027 

(0.079) 

19.091 

(0.061) 

18.932 

(0.054) 

18.941 

(0.073) 

18.9121 

For  both  the  AR(1)  and  EAR(l)  processes,  the  covariance  function  Rk  =  k  = 
0,  ±1,  ±2,  —  The  covariance  function  of  the  M/M/1  waiting  time  process  is  more  com¬ 
plicated  (cf.  Daley  [5]). 

We  simulated  the  above  stochastic  processes  over  a  variety  of  parameter  values;  rep¬ 
resentative  results  are  presented  in  Table  I  (AR{1)  with  ^  =  0.9),  Table  2  (EAR(l)  with 
^  =  0.9),  and  Table  3  (M/M/1  waiting  time  process  with  arrival  rate  0.8  and  service  rate 
1.0).  Each  table  entry  in  a  row  is  based  on  the  same  100,000  independent  replications 
of  the  stochastic  process.  The  number  in  parentheses  below  an  entry  is  the  standard 
error  of  that  entry.  Each  of  the  replications  was  initialized  from  the  appropriate  steady- 
state  distribution.  All  uniform  [normal]  random  variates  were  generated  from  algorithm 
UNIF  [TRPNRM]  in  Bratley,  Fox,  and  Schrage  [3];  exponential  deviates  used  inversion;  the 
M/M/1  waiting-time  process  was  generated  from  an  algorithm  due  to  Schmeiser  [18]. 
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Table  2;  Estimated  Expected  Values  of  Various  Variance  Estimators  —  EAR(l),  4>  =  0.9. 
(Note  that  <t*  =  19.0  for  this  process.) 


1  n 

A^ifr.n) 

C(go',  n) 

C(gAD‘,n) 

C(gV,n) 

C(ff;:n) 

C(?6;»») 

H 

HiBlI 

iiftil 

IHI 

0.322 

(0.005) 

0.299 

(0.004) 

0.226 

(0.003) 

0.303 

(0.004) 

3.5245 

8 

IHI 

im 

■HI 

0.946 

(0.010) 

0.764 

(0.007) 

0.677 

(0.006) 

0.897 

(0.009) 

6.1855 

16 

2.859 

(0.021) 

IHI 

HI 

1.717 

(0.011) 

2.770 

(0.021) 

2.188 

(0.015) 

1.866 

(0.012) 

2.627 

(0.019) 

9.8346 

32 

liil 

4.160 

(0.022) 

Biil 

Bin 

13.5681 

64 

HUM 

11.934 

(0.062) 

9.932 

(0.047) 

8.570 

(0.039) 

11.388 

(0.058) 

16.1908 

128 

BSiH 

13.218 

(0.053) 

11.536 

(0.044) 

16.198 

(0.077) 

14.491 

(0.059) 

12.884 

(0.050) 

15.602 

(0.072) 

17.5938 

256 

fl^M 

IH 

15.821 

(0.055) 

14.289 

(0.047) 

18.081 

(0.080) 

17.246 

(0.063) 

16.138 

(0.055) 

17.629 

(0.074) 

18.2969 

512 

17.345 

(0.056) 

16.161 

(0.048) 

18.771 

(0.080) 

18.511 

(0.064) 

17.925 

(0.057) 

18.486 

(0.074) 

18.6484 

1024 

18.374 

(0.084) 

18.841 

(0.086) 

mijijiji 

■PfihylM 

18.8242 

2048 

Ri9 

18.9121 

i 
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Table  3:  Estimated  Expected  Values  of  Various  Variance  Estimators  —  M/M/1  Waiting 
Time  Process;  arrival  rate  =  0.8,  service  rate  =  1.0 


C{g 

21.6 

18.0 

29.1 

22.9 

19.5 

27.6 

(0.1) 

(0.1) 

(0.2) 

(0.1) 

(0.1) 

(0-2) 

66.4 

55.5 

89.5 

70.5 

59.7 

85.1 

(0.4) 

(0.3) 

(0.5) 

(0.4) 

(0.3) 

(0.5) 

1940.9 

1870.6 

2008.2 

2009.8  ] 

(8-9) 

(8.5) 

(11.1) 

(10.3)  1 

6  Discussion 


This  section  summarizes  and  discusses  the  exact  and  estimated  expectation  and  variance 
results  for  the  variance  estimators  examined  in  §§2-5.  Recall  that  we  obtained  exact 
results  for  area  estimators  in  §2  and  for  CvM  estimators  in  §3.  We  also  gave  exact  results 
for  a  specific  stochastic  process,  the  MA(1),  in  §4  and  empirical  results  for  the  AR(1), 
EAR(l),  and  M/M/1  in  §5.  The  estimated  expected  values  given  in  Tables  1-3  are  based 
on  100,000  replications;  thus,  one  can  obtain  estimated  variance  results  from  the  Monte 
Carlo  runs  by  squaring  the  estimated  standard  errors  (in  parentheses  below  the  estimated 
expected  values)  and  then  multiplying  by  100,000. 

For  each  of  the  stochastic  processes  under  study,  the  expected  value  of  the  unweighted 
area  estimator  A*(/o;n)  converged  relatively  slowly  to  as  n  increased.  This  phe¬ 
nomenon  is  due  to  the  estimator’s  comparatively  high  Bias(A*(/o;  n))  %  O-y/n  (Example 
1;  also  see  [17]). 

The  expected  value  of  the  unweighted  CvM  estimator  C(po; «)  converged  more  slowly 
than  that  of  A*(/o;n)  to  This  makes  sense  since  Bias(C(i?o; «))  ^  S^/n  (Example  3) 
is  higher  than  A^(/o;n)’s.  On  the  plus  side,  we  see  that  for  large  n, 

0  Ra* 

Var(C(go;n))/Var(A^(/o;n))  «  =  0.4, 

as  suggested  by  Theorem  1  and  Example  6. 

The  expected  value  of  the  Anderson-Darling  CvM  estimator  C(gAD;  n)  converged  even 
more  slowly  than  that  of  C(goin)  to  o’.  Although  we  did  not  prove  the  inferiority  of 
the  convergence  rate  of  E[C(^/4o;n)|  to  o’  for  general  stationary  processes,  the  evidence 
provided  by  Example  14  and  the  empirical  work  seems  to  be  overwhelming.  Some  calculus 
shows  that 

Var(C(gAD))  =  (j  -  «  0.5797o^ 

(cf.  [21]);  but  it  is  of  little  comfort  that  C(gAD',  n)  has  the  lowest  variance  of  the  estimators 
under  study. 

The  expected  values  of  the  first-order  unbiased  quadratic  area  estimator  for  o’, 
A’(/2;n),  and  the  minimum-variance  first-order  unbiased  quadratic  CvM  estimator  for 
o’,  C(g2i  n),  converged  comparatively  quickly  to  o’  as  n  increased;  the  rapid  convergence 
is  a  direct  consequence  of  the  first-order  unbiasedness  of  the  estimators.  For  large  n,  we 
see  from  Examples  12  and  13  and  the  empirical  tables  that 

Var(C(^J;n))/Var(A’(/2;n))  as  121/140, 

as  predicted  by  Theorem  1  and  Example  7. 

The  minimum-variance  first-order  unbiased  fourth-  and  sixth-degree  CvM  estimators, 
respectively,  possess  expected  values  that  converge  to  o’  almost 
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(but  not  quite)  as  quickly  as  those  of  A^{f2\ n)  and  C{gy,  n).  A  favorable  property  of  these 
higher-degree  estimators  is  that  they  have  reduced  variances;  the  variance  improvements 
are  along  the  lines  described  in  Example  9. 

Recall  that  the  estimator  C(gj;n)  is  first-order  unbiased  for  (Example  5).  We 
can  compare  the  estimated  expected  values  from  Tables  1  and  2  with  the  corresponding 
actual  <7|[-vaiues  (given  in  the  last  column  of  the  tables).  We  see  that  the  bias  of  C(9|;  n) 
as  an  estimator  for  a\  is  about  the  same  as  Bias(/l^(/2;n))  and  Bias(C(y2;n));  fhis  is 
particularly  true  for  large  n.  In  addition,  Var(C(p2;*i))  only  a  little  smaller  than 
Var(C(</2;”))-  Thus,  we  do  not  seem  to  gain  much  by  using  C(55;n)  to  estimate 

The  bottom  line:  Of  the  estimators  studied  so  far,  it  appears  that  C{g2'.n)  performs  the 
best. 


7  Conclusions 

In  this  article,  we  introduced  a  class  of  estimators  for  —  limn-.oo  nVar(F„)  that  are 
similar  to  Cramer-von  Mises  statistics.  Using  appropriate  weighting  functions,  our  CvM 
estimators  were  shown  to  be  first-order  unbiased,  and  asymptotic  variance  reductions 
of  up  to  60%  (compared  to  the  weighted  area  estimator)  were  achievable.  Further,  the 
variance-optimal  weighting  functions  can  be  computed  independently  of  the  output  pro¬ 
cess.  Our  conclusions  were  supported  with  an  analytical  example  using  the  MA(1)  pro¬ 
cess. 

Although  the  estimators  are  all  asymptotically  unbiased  for  finite-sample  bias 
can  be  substantial.  Analytical  and  empirical  work  showed  that  the  bias  of  the  CvM 
estimators  converged  to  zero  at  least  as  fast  as  that  of  the  weighted  area  estimators,  and 
the  CvM  estimators  had  smaller  variance.  Thus,  if  the  sample  size  is  sufficiently  large, 
the  CvM  estimators  proved  to  be  more  efficient  than  the  weighted  area  estimators. 

As  discussed  in  [10],  it  is  possible  to  augment  the  basic  CvM  variance  estimator  in  a 
number  of  ways. 

1.  One  csui  show  that  the  unweighted  CvM  and  area  estimators  are  highly  correlated. 
This  suggests  that  certain  linear  combinations  of  the  area  and  CvM  estimators  will 
give  rise  to  estimators  having  comparatively  lower  variance. 

2.  All  of  our  work  so  far  has  assumed  that  we  have  one  long  batch  of  n  observations. 
An  alternative  way  of  organizing  the  data  is  to  break  the  n  observations  into  6 
contiguous,  nonoverlapping  batches,  each  of  size  m  (assume  n  =  frm).  This  leads  to 
another  interesting  problem — that  of  examining  the  consequences  of  batching  the 
data  and  then  forming  CvM  estimators  from  each  batch.  Intuitively,  batching  of  the 
data  will  tend  to  increase  estimator  bias  while  decreasing  estimator  variance — of 
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course,  one  can  quantify  the  trade-off  by  calculating  the  batched  CvM  estimator’s 
mean  '■quared  error.  (See  Schmeiser  and  Song  [19].) 

3.  We  can  also  apply  the  methodology  of  Meketon  and  Schmeiser  [15]  in  which  the 
n  observations  are  broken  into  n  —  m  +  1  overlapping  batches,  each  of  size  m. 
Although  the  bias  of  the  resulting  overlapped  CvM  estimator  is  the  same  as  that  of 
the  batched  CvM  estimator,  the  overlapped  estimator’s  variance  is  much  smaller. 

The  above  problems  are  the  subjects  of  ongoing  research  and  will  be  the  topics  of  a 
future  paper. 
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Appendix 

This  appendix  contains  the  proof  of  Theorem  2.  Before  proving  the  theorem,  we  state 
and  prove  a  series  of  lemmas.  First,  we  define  the  cumulative  sums  Z*  s  JZjsi  Y> 
the  variance  time  curve  ^(1:)  =  Var(Zfc),  k  =  1,2,  ...,n  (see  [11]).  Since  g{t)  is  assumed 
to  be  continuous  and  bounded  on  [0, 1],  we  denote  M  =  supo<Ki  1^(01  <  oo. 

Lemma  2  Under  the  Assumptions  of  §3, 

V{n)  =  n<7* -f  7  -  2  ^(n  -  t)i2,  =  tkt*  +  7 -f  o(l). 

isn 

Proof  Follows  from  the  arguments  leading  to  Equation  (2).  □ 


Lemma  3  (The  discrete-time  version  of  a  result  given  in  [11].)  Under  the  Assumptions 
of  §3,  if  n  >  ib,  we  have 

2  Cov(Zn,  Z*)  =  V(n)  -h  V{k)  -  V{n  -  it). 

Proof  By  stationary  increments, 

V(n)  =  Var(Z(n)  -  Z(ib) -b  Z(ib)) 

=  Var(Z(n)  -  Z{k))  -b  Var(Z(ib))  +  2Cov(Z(n)  -  Z(/b),  Z{k)) 

=  Var(Z(n-ib))-b  V(ib)-b2Cov(Z(n),Z(ib))-2V(jb).  □ 
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Lemma  4  Under  the  Assumptions  of  §3, 


=  o(n’). 


Lemma  5  Under  the  Assumptions  of  §3, 


-  k)Rj-  U-in-k))Rj 

*=1  jafc  jsn-fc 

<  jw  i:  1 13  f  +  “f; 

kxl  jsn  jsk  jsn-k 


L  j="  *»i 

=  o(n^). 


Lemma  6  Under  the  Assumptions  of  §3, 

k»\  ”  jsk 


k=l jsk 

=  +  f; 

**l  >*n+I 

=  o(n). 


Lemma  7  Under  the  Assumptions  of  §3, 

+ 7 +  <Kl))L  «(“)*’  *  («<^’  +  7)Sl9(^)*’  +  o(n^). 


Proof  Follows  from  Lemmas  2  and  4.  □ 


20 


Lemma  8  Under  the  Assumptions  of  §3, 

f;<,(i)*cov(z„2.) 

=  Y.9(-)I‘  *0’ +  I  +  Di -»)«,  + Dj -*)«,  -  E  (i  -  (n  -  *))«, 

fcal  ^  J=n-k 

=  + +  <>(”*)• 
ksl  ”  *=1  ” 

Proof  Follows  from  Lemmas  2,  3,  and  5.  □ 

Lemma  9  Under  the  Assumptions  of  §3, 

k=l  ”  *=l  ”  [  ja* 

k^x  ^  t^x  ^ 

Proof  Follows  from  Lemmas  2  and  6.  O 


Lemma  10  (Trapezoid  Rule).  If  e''(f)  is  continuous  and  bounded  V<  €  [0, 1],  then 


/*  /  \  j  1  ^  /*\  «(0)  “■«(!)  . 


o{l/n). 


We  can  finally  prove  Theorem  2. 
Proof  (of  Theorem  2). 


=  ^  -  h  p(^‘Cc(z„z.) + i  |:,(i)v(*) 

(by  Lemmas  7,  8,  and  9,  and  algebra) 

=  1 9m  -t^)dt  +  ^  ['  gm^  -«+!)*  +  o(l/n) 

Jo  n  Jo 


(by  Lemma  10). 


Application  of  Assumption  6  completes  the  proof.  □ 


Remark  3  It  is  quite  a  bit  easier  to  prove  the  continuous-time  version  of  the  theorem 

(cf.  (111). 
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